subroutine  compute_Hutton_Boore_ML(c,nd,stanum,evnum,nsta,nev,sigma2_HB)
    implicit none
    
    include 'seisev_EVENT_PHASE_structure.FI'
    include 'inversione_pesata_multipar.FI'
    
    integer i,k,nd,isw,nv
    integer stanum(150000),evnum(150000),nsta,nev
    integer first,last,dst
    real*4  dif,r
    real*4  v(5000),wmag,wsd,wmax,wmin,ww(5000),ruse
    real*8  sigma2_HB
    data isw/0/

! compute all the station Hutton-Boore magnitudes    
    do i=1,nd       
      r = c(i).s_h_hypdist
      ! compute station magnitude
      c(i).MLHB_sta = log10(c(i).WA_pha_amp) + 1.110 *(r /100.) +  &
                      0.00189*(r - 100.) + 3.0
    enddo
    
!   Prima ricalcolo le magnitudo di evento con la mediana,
!   poi con la media di Huber (tanto per vedere la differenza
!   delle due varianze).
    
    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1
    
    sigma2_HB=0.0d00
    do while (last .le. nd  .and. first .le. nd)
    
      k= 0
      do i=first,last
          k=k+1 
          v(k) = c(i).MLHB_sta
      enddo
    
      call mediana(v,k,wmag)      
    
      k=0
      do i=first,last
          k=k+1
          c(i).MLHB_ev = wmag     !event mag
          dif = c(i).MLHB_sta - wmag
          sigma2_HB = sigma2_HB + dif*dif
      enddo   
                 
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1
    
    enddo        
    
    write(*,*)' calcolo delle ML Hutton-Boore con la mediana. Sigma2 = ',sigma2_HB    
    
    ! go throught the events and compute the event HB Magnitudes.
    ! ora calcolo le magnitudo HB finali con Huber

    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1

    sigma2_HB=0.0d00
    do while (last .le. nd  .and. first .le. nd)

      k= 0
      do i=first,last
          k=k+1 
          v(k) = c(i).MLHB_sta
      enddo

      ruse = 0.3
      wmag = 0.
      wsd = 0.
      call Huber_dev(k,v,ruse,wmag,wsd,wmax,wmin,ww)
      ! memorizza i risultati   !
      !
      
      k=0
      do i=first,last
          k=k+1
          c(i).MLHB_ev = wmag     !event mag
          dif = c(i).MLHB_sta - wmag
          sigma2_HB = sigma2_HB + dif*dif
      enddo   
                 
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1

    enddo        


    write(*,*)' calcolo delle ML HB con Huber. Sigma2 = ',sigma2_HB    


return
end
    
    
 